Code finally compiles as of 414eac6!

full_integration passes as of 83f5168, which does not test the fct at all, but indicates that the new li_mesh and li_config modules are operational.

Testing on a364f7a shows that the code will run with config_thickness_advection = 'fo' and config_tracer_advection = 'fct', but from examining damage output, tracers don't seem to actually be advecting.

As of e85d3fd, fct tracer advection is at least doing something, although is not yet verified. Here are three MISMIP runs after 100 years with velo solver turned off, no sub-shelf melt, calculate-damage = true. All use 'fo' thickness advection. Left: no tracer advection; Middle: fo tracer advection; right: fct tracer advection:

White contour is damage = 0.4 from the 'fo' case.

As of 8734303, fct thickness advection is operational, and seems to be working correctly. One issue is that it seems to require a very small CFL fraction to avoid terracing effects.
From left to right: fo thickness advection with CFL fraction = 0.8; fct thickness advection w/ CFL fraction 0.8; fct thickness advection w/ CFL fraction 0.1:

CFL fraction of 0.1 (~40 day time steps) seems to be sufficiently low here. I tested another run w/ a CFL fraction of 0.01 (~4 days) and the differences at year 100 are generally <1 m:

What we are seeings seems to be a known issue w/ FCT: "Another infamous byproduct of FCT manifests itself in distortions of a smooth profile. This phenomenon is known as terracing and represents ‘an integrated, nonlinear effect of residual phase errors’ [64] or, loosely speaking, ‘the ghosts of departed ripples’ [7]. A particularly severe form of terracing is caused by the linear instability of the high-order scheme. For this reason, we do not recommend the use of the forward Euler time-stepping (θ = 0) even though the flux-corrected scheme proves positivity-preserving under the CFL-like time step restriction (62)."
(found on page 178 of https://link.springer.com/content/pdf/10.1007/978-94-007-4038-9.pdf)

Adding to comment above: CFL fraction of 0.4 also drastically reduces the terracing and gives a thickness field visually indistinguishable from CFL fractions of 0.1 and 0.01, so that's great news! Differences between runs w/ CFL fractions of 0.4 and 0.1:

Just for comparison, here are the thickness differences at year 100 when using FO thickness advection with CFL of 0.8 vs 0.1:

So the FCT solution is definitely more sensitive to the time step than the FO solution, even when the time step is quite small.

Testing on Thwaites mesh indicates issue with FCT thickness advection during advance, which is not surprising. It turns out that boundaryCell is not being calculated, and thus neither is advMaskHighOrder. This leads to very high velocities and strange geometries in areas where ice advances beyond the initial extent.

I have confirmed that results are BFB for restarts and decomposition tests (64 and 128 cores) using the slotted cylinder test case. I should also put together a version of full_integration that uses fct.
Here are the restart test results, with the full run on top, a restart run branching from year 10 in the middle and the difference at year 25 on bottom. This uses 3rd order, beta = 0.5, and CFL fraction = 0.25 (hence the ugliness)

And here is the same for the decomp test. Top is 128 procs, middle is 64 procs, bottom is difference:

Here are results at year 25 with CFL fractions of 0.25 (top), 0.1 (middle) and 0.05 (bottom):

Uh oh, it appears that passiveTracer is not being conserved in test runs with fct. I ran the notched slab test that I was using above for 100 years and summed up passiveTracer at each time step: ncap2 -s 'tracer_sum = (passiveTracer * thickness * areaCell).sum($nCells)' output_all_timesteps.nc tmp.nc.
fct:

fo conserves just fine:

Both runs do conserve mass, however. Blue is fo, orange is fct:

Update:
I rewound the branch to 4c87809 and the issue is still there, but it is less egregious and drifts the opposite way:

8d2f023 fixed the tracer conservation issue detailed above! It still doesn't conserve to machine precision, but a 100 year test run on the grounded slab test case showed about 0.06% drift by the end of the simulation. That seems pretty close!

Compare this with the drift above

@dengwirda, see note about slight non-conservation of tracers in Trevor's current implementation (#70 (comment)). Should we be expecting machine precision conservation? If so, do you think we should worry that we're not achieving that given how small the non-conservation is?

I suppose we would be introducing energy conservation errors, given that we will be advecting temperature with this scheme.

Testing

@trhille , thanks for making the final (more complicated than expected) adjustment. I tested one more time and all tests pass:

Test Runtimes:
00:19 PASS landice_dome_2000m_sia_restart_test
00:19 PASS landice_dome_2000m_sia_decomposition_test
00:16 PASS landice_dome_variable_resolution_sia_restart_test
00:06 PASS landice_dome_variable_resolution_sia_decomposition_test
00:26 PASS landice_enthalpy_benchmark_A
00:24 PASS landice_eismint2_decomposition_test
00:23 PASS landice_eismint2_enthalpy_decomposition_test
00:29 PASS landice_eismint2_restart_test
00:28 PASS landice_eismint2_enthalpy_restart_test
00:22 PASS landice_greenland_sia_restart_test
00:12 PASS landice_greenland_sia_decomposition_test
00:20 PASS landice_hydro_radial_restart_test
00:12 PASS landice_hydro_radial_decomposition_test
00:44 PASS landice_humboldt_mesh-3km_decomposition_test_velo-none_calving-none_subglacialhydro
00:30 PASS landice_humboldt_mesh-3km_restart_test_velo-none_calving-none_subglacialhydro
00:26 PASS landice_dome_2000m_fo_decomposition_test
00:30 PASS landice_dome_2000m_fo_restart_test
00:11 PASS landice_dome_variable_resolution_fo_decomposition_test
00:16 PASS landice_dome_variable_resolution_fo_restart_test
00:18 PASS landice_circular_shelf_decomposition_test
00:44 PASS landice_greenland_fo_decomposition_test
00:44 PASS landice_greenland_fo_restart_test
00:41 PASS landice_thwaites_fo_decomposition_test
00:32 PASS landice_thwaites_fo_restart_test
00:15 PASS landice_thwaites_fo-depthInt_decomposition_test
00:22 PASS landice_thwaites_fo-depthInt_restart_test
00:58 PASS landice_humboldt_mesh-3km_restart_test_velo-fo_calving-von_mises_stress_damage-threshold_faceMelting
00:32 PASS landice_humboldt_mesh-3km_restart_test_velo-fo-depthInt_calving-von_mises_stress_damage-threshold_faceMelting
Total runtime 12:04
PASS: All passed successfully!

I will also test the FCT tests you have in MPAS-Dev/compass#650

Testing, Part 2

The FCT suite all passes. I did not look at any of the output, but they all run and pass validation. For the sake of merging this PR as a checkpoint, I think that is sufficient.

Test Runtimes:
00:19 PASS landice_dome_2000m_fo_fct_decomposition_test
00:35 PASS landice_dome_2000m_fo_fct_restart_test
00:16 PASS landice_dome_variable_resolution_fo_fct_decomposition_test
00:14 PASS landice_dome_variable_resolution_fo_fct_restart_test
00:42 PASS landice_greenland_fo_fct_decomposition_test
00:47 PASS landice_greenland_fo_fct_restart_test
00:38 PASS landice_thwaites_fct_decomposition_test
00:34 PASS landice_thwaites_fct_restart_test
00:35 PASS landice_humboldt_mesh-3km_restart_test_velo-fo_advec-fct_calving-von_mises_stress_damage-threshold_faceMelting
Total runtime 04:46
PASS: All passed successfully!